module stepon
!----------------------------------------------------------------------- 
! 
! Purpose: 
! Module for time-stepping of the CAM Eulerian Spectral dynamics.
! 
!-----------------------------------------------------------------------
  use shr_kind_mod,     only: r8 => shr_kind_r8
  use shr_sys_mod,      only: shr_sys_flush
  use pmgrid,           only: plev, plat, plevp, plon, beglat, endlat
  use spmd_utils,       only: masterproc
  use scanslt,          only: advection_state
  use prognostics,      only: ps, u3, v3, t3, q3, qminus, div, &
                              dpsl, dpsm, omga, phis, n3, n3m2, n3m1
  use camsrfexch_types, only: cam_out_t     
  use ppgrid,           only: begchunk, endchunk
  use physics_types,    only: physics_state, physics_tend
  use time_manager,     only: is_first_step
  use infnan,           only: nan
  use iop,              only: setiopupdate, readiopdata
  use scamMod,          only: use_iop,doiopupdate,use_pert_frc,wfld,wfldh,single_column
  use perf_mod
  implicit none

  private   ! By default make all data and methods private to this module
!
! Public methods
!
  public stepon_init     ! Initialization
  public stepon_run1     ! Run method phase 1
  public stepon_run2     ! Run method phase 2
  public stepon_run3     ! Run method phase 3
  public stepon_final    ! Finalization
!
! Private module data
!
  save
  type(physics_state), pointer :: phys_state(:)   ! Physics state data
  type(physics_tend ), pointer :: phys_tend(:)    ! Physics tendency data

  real(r8) :: detam(plev)               ! intervals between vert full levs.
  real(r8) :: cwava(plat)               ! weight applied to global integrals
  real(r8), allocatable :: t2(:,:,:)    ! temp tendency
  real(r8), allocatable :: fu(:,:,:)    ! u wind tendency
  real(r8), allocatable :: fv(:,:,:)    ! v wind tendency
  real(r8), allocatable :: flx_net(:,:) ! net flux from physics
  real(r8), allocatable :: fq(:,:,:,:)  ! Q tendencies,for eul_nsplit>1
  real(r8), allocatable :: t2_save(:,:,:)    ! temp tendency
  real(r8), allocatable :: fu_save(:,:,:)    ! u wind tendency
  real(r8), allocatable :: fv_save(:,:,:)    ! v wind tendency
  real(r8) :: coslat(plon)              ! cosine of latitude
  real(r8) :: rcoslat(plon)             ! Inverse of coseine of latitude
  real(r8) :: rpmid(plon,plev)          ! inverse of midpoint pressure
  real(r8) :: pdel(plon,plev)           ! Pressure depth of layer
  real(r8) :: pint(plon,plevp)          ! Pressure at interfaces
  real(r8) :: pmid(plon,plev)           ! Pressure at midpoint
  real(r8) :: dtime = nan               ! timestep size
  type(advection_state) :: adv_state    ! Advection state data

!======================================================================= 
contains
!======================================================================= 

!
!======================================================================= 
!

subroutine stepon_init( gw, etamid, dyn_in, dyn_out )
!----------------------------------------------------------------------- 
! 
! Purpose:  Initialization, primarily of dynamics.
!
!----------------------------------------------------------------------- 
   use dyn_comp,       only: dyn_import_t, dyn_export_t
   use scanslt,        only: scanslt_initial
   use commap,         only: clat
   use constituents,   only: pcnst
   use physconst,      only: gravit
   use rgrid,          only: nlon
   use hycoef,         only: hyam, hybm
   use time_manager,   only: get_step_size
   use eul_control_mod,only: eul_nsplit
#if ( defined BFB_CAM_SCAM_IOP )
   use iop,            only:init_iop_fields
#endif
!-----------------------------------------------------------------------
! Arguments
!
  real(r8), intent(out) :: gw(plat)                  ! Gaussian weights
  real(r8), intent(out) :: etamid(plev)              ! vertical coords at midpoints
  type(dyn_import_t) :: dyn_in                       ! included for compatibility
  type(dyn_export_t) :: dyn_out                      ! included for compatibility
!-----------------------------------------------------------------------
!  Local variables
!
   integer :: k, lat, i
   !-----------------------------------------------------------------------

   call t_startf ('stepon_startup')

   dtime = get_step_size()
   !
   ! Define eta coordinates: Used for calculation etadot vertical velocity 
   ! for slt.
   !
   do k=1,plev
      etamid(k) = hyam(k) + hybm(k)
   end do

   call scanslt_initial( adv_state, etamid, gravit, gw, detam, cwava )
   !
   ! Initial guess for trajectory midpoints in spherical coords.
   ! nstep = 0:  use arrival points as initial guess for trajectory midpoints.
   ! nstep > 0:  use calculated trajectory midpoints from previous time 
   ! step as first guess.
   ! NOTE:  reduce number of iters necessary for convergence after nstep = 1.
   !
   if (is_first_step()) then
      do lat=beglat,endlat
	 if (.not. single_column) then
            do i=1,nlon(lat)
               coslat(i) = cos(clat(lat))
               rcoslat(i) = 1._r8/coslat(i)
            end do
         endif
         !     
         ! Set current time pressure arrays for model levels etc.
         !
         call plevs0(nlon(lat), plon, plev, ps(1,lat,n3), pint, pmid, pdel)
         !
         do k=1,plev
            do i=1,nlon(lat)
               rpmid(i,k) = 1._r8/pmid(i,k)
            end do
         end do

	 if (.not. single_column) then
         !
         ! Calculate vertical motion field
         !
            call omcalc (rcoslat, div(1,1,lat,n3), u3(1,1,lat,n3), v3(1,1,lat,n3), dpsl(1,lat), &
                      dpsm(1,lat), pmid, pdel, rpmid   ,pint(1,plevp), &
                      omga(1,1,lat), nlon(lat))
         else
         
            omga(1,:,lat)=wfld(:)
         endif
      end do
   end if

   allocate(t2(plon,plev,beglat:endlat))
   allocate(fu(plon,plev,beglat:endlat))
   allocate(fv(plon,plev,beglat:endlat))
   allocate( flx_net(plon,beglat:endlat))
   if (eul_nsplit>1) then
      allocate(fq(plon,plev,pcnst,beglat:endlat))
      allocate(t2_save(plon,plev,beglat:endlat))
      allocate(fu_save(plon,plev,beglat:endlat))
      allocate(fv_save(plon,plev,beglat:endlat))
   endif
   !
   ! Beginning of basic time step loop
   !
   call t_stopf ('stepon_startup')


#if ( defined BFB_CAM_SCAM_IOP )
   if (is_first_step()) then
      call init_iop_fields()
   endif
#endif
end subroutine stepon_init

!
!======================================================================= 
!

subroutine stepon_run1( ztodt, phys_state, phys_tend , dyn_in, dyn_out)
!----------------------------------------------------------------------- 
! 
! Purpose:  Phase 1 run method of dynamics. Set the time-step
!           to use for physics. And couple from dynamics to physics.
!
!----------------------------------------------------------------------- 
  use dyn_comp,       only: dyn_import_t, dyn_export_t
  use time_manager,   only: get_nstep
  use prognostics,    only: pdeld
  use phys_buffer,    only: pbuf
  use dp_coupling,    only: d_p_coupling
  use eul_control_mod,only: eul_nsplit

  real(r8), intent(out) :: ztodt            ! twice time step unless nstep=0
  type(physics_state), intent(inout) :: phys_state(begchunk:endchunk)
  type(physics_tend),  intent(out)   :: phys_tend(begchunk:endchunk)
  type(dyn_import_t) :: dyn_in                       ! included for compatibility
  type(dyn_export_t) :: dyn_out                      ! included for compatibility

  ztodt = 2.0_r8*dtime
  !
  ! If initial time step adjust dt
  !
  if (is_first_step()) ztodt = dtime

  ! subcycling case, physics dt is always dtime
  if (eul_nsplit>1) ztodt = dtime	

  !
  ! adjust hydrostatic matrices if the time step has changed.  This only
  ! happens on transition from time 0 to time 1. 
  ! (this is now initialized in dynpgk)
  !if (get_nstep() == 1) then
  !   call settau(dtime)
  !end if
  !
  ! Dump state variables to IC file
  !
  call t_startf ('diag_dynvar_ic')
  call diag_dynvar_ic (phis, ps(1,beglat,n3m1), t3(1,1,beglat,n3m1), u3(1,1,beglat,n3m1), &
                       v3(1,1,beglat,n3m1), q3(1,1,1,beglat,n3m1) )
  call t_stopf ('diag_dynvar_ic')
  !
  !----------------------------------------------------------
  ! Couple from dynamics to physics
  !----------------------------------------------------------
  !
  call t_startf ('d_p_coupling')
  call d_p_coupling (ps(1,beglat,n3m2), t3(1,1,beglat,n3m2), u3(1,1,beglat,n3m2), &
                     v3(1,1,beglat,n3m2), q3(1,1,1,beglat,n3m2), &
                     omga, phis, phys_state, phys_tend, pbuf, pdeld(:,:,:,n3m2))
  call t_stopf  ('d_p_coupling')
end subroutine stepon_run1

!
!======================================================================= 
!

subroutine stepon_run2( phys_state, phys_tend, dyn_in, dyn_out )
!----------------------------------------------------------------------- 
! 
! Purpose:  Phase 2 run method of dynamics. Couple from physics
!           to dynamics.
!
!----------------------------------------------------------------------- 
  use dyn_comp,       only: dyn_import_t, dyn_export_t
  use dp_coupling,    only: p_d_coupling
  type(physics_state), intent(in):: phys_state(begchunk:endchunk)
  type(physics_tend), intent(in):: phys_tend(begchunk:endchunk)
  type(dyn_import_t) :: dyn_in                       ! included for compatibility
  type(dyn_export_t) :: dyn_out                      ! included for compatibility

  call t_startf ('p_d_coupling')
  call p_d_coupling (phys_state, phys_tend, t2, fu, fv, flx_net, &
                     qminus(1,1,1,beglat) )
  call t_stopf  ('p_d_coupling')
end subroutine stepon_run2

!
!======================================================================= 
!

subroutine stepon_run3( ztodt, etamid, cam_out, phys_state, dyn_in, dyn_out )
!----------------------------------------------------------------------- 
! 
! Purpose:  Final phase of dynamics run method. Run the actual dynamics.
!
!----------------------------------------------------------------------- 
  use dyn_comp,       only: dyn_import_t, dyn_export_t
  use eul_control_mod,only: eul_nsplit
  real(r8), intent(in) :: ztodt            ! twice time step unless nstep=0
  type(cam_out_t), intent(inout) :: cam_out(begchunk:endchunk)
  real(r8), intent(in) :: etamid(plev)     ! vertical coords at midpoints
  type(physics_state), intent(in):: phys_state(begchunk:endchunk)
  type(dyn_import_t) :: dyn_in                       ! included for compatibility
  type(dyn_export_t) :: dyn_out                      ! included for compatibility
  real(r8) :: dt_dyn0,dt_dyn 
  integer :: stage
  if (single_column) then
     
     ! Determine whether it is time for an IOP update;
     ! doiopupdate set to true if model time step > next available IOP
     if (use_iop) then
        call setiopupdate
     end if
     
     ! Update IOP properties e.g. omega, divT, divQ
     
     if (doiopupdate) call readiopdata()
     
  endif

  !----------------------------------------------------------
  ! DYNPKG Call the Dynamics Package
  !----------------------------------------------------------
  call t_startf ('dynpkg')

  if (eul_nsplit==1) then	
     call dynpkg(adv_state, t2      ,fu      ,fv      ,etamid  ,  &
       cwava   ,detam   ,flx_net ,ztodt)
  else
     dt_dyn0 = ztodt/eul_nsplit
     dt_dyn = dt_dyn0
     if (is_first_step()) dt_dyn = 2*dt_dyn0

     ! convert q adjustment to a tendency
     fq = (qminus(:,:,:,:) - q3(:,:,:,:,n3m2))/ztodt
     ! save a copy of t2,fu,fv
     t2_save=t2		
     fu_save=fu
     fv_save=fv

     call apply_fq(qminus,q3(:,:,:,:,n3m2),fq,dt_dyn0)
     call dynpkg(adv_state, t2      ,fu      ,fv      ,etamid  ,  &
       cwava   ,detam   ,flx_net ,dt_dyn0)

     do stage=2,eul_nsplit	
        t2=t2_save		
        fu=fu_save
        fv=fv_save
        call apply_fq(qminus,q3(:,:,:,:,n3m2),fq,dt_dyn)
        call dynpkg(adv_state, t2      ,fu      ,fv      ,etamid  ,  &
             cwava   ,detam   ,flx_net ,dt_dyn)
     enddo
  endif

  call t_stopf  ('dynpkg')
end subroutine stepon_run3



subroutine apply_fq(qminus,q3,fq,dt)
   use shr_kind_mod, only: r8 => shr_kind_r8
   use pmgrid,       only: plon, plat, plev, plevp, beglat, endlat
   use rgrid,        only: nlon
   use constituents,   only: pcnst

   real(r8), intent(in) :: q3(plon,plev,beglat:endlat,pcnst)
   real(r8), intent(in) :: fq(plon,plev,beglat:endlat,pcnst)
   real(r8), intent(out) :: qminus(plon,plev,beglat:endlat,pcnst)
   real(r8), intent(in) :: dt

   !local 	
   real(r8) :: q_tmp,fq_tmp
   integer :: q,c,k,i

   do q=1,pcnst 
   do c=beglat,endlat
   do k=1,plev
   do i=1,nlon(c)
      fq_tmp = dt*fq(i,k,c,q)
      q_tmp  = q3(i,k,c,q)
      ! if forcing is > 0, do nothing (it makes q less negative)
      if (fq_tmp<0 .and. q_tmp+fq_tmp<0 ) then
         ! reduce magnitude of forcing so it wont drive q negative 
         ! but we only reduce the magnitude of the forcing, dont increase
         ! its magnitude or change the sign
         
         ! if q<=0, then this will set fq=0  (q already negative)
         ! if q>0, then we know from above that fq < -q < 0, so we 
         ! can reduce the magnitive of fq by setting fq = -q:
         fq_tmp = min(-q_tmp,0d0)
      endif
      qminus(i,k,c,q) = q_tmp + fq_tmp
   enddo
   enddo
   enddo
   enddo
   	
end subroutine


!
!======================================================================= 
!

subroutine stepon_final(dyn_in, dyn_out)
!----------------------------------------------------------------------- 
! 
! Purpose:  Stepon finalization.
!
!----------------------------------------------------------------------- 
   use dyn_comp,       only: dyn_import_t, dyn_export_t
   use scanslt, only: scanslt_final
   type(dyn_import_t) :: dyn_in                       ! included for compatibility
   type(dyn_export_t) :: dyn_out                      ! included for compatibility

   call print_memusage ('End stepon')
   call scanslt_final( adv_state )
   deallocate(t2)
   deallocate(fu)
   deallocate(fv)
   deallocate(flx_net)

end subroutine stepon_final
!
!======================================================================= 
!

end module stepon
